Demonstration: Reading Data¶
This notebook explains the data output from the GEODYN Starlette-SLR run.
First, load the data using the reading functions.
[1]:
%load_ext autoreload
%autoreload 2
import pandas as pd
import numpy as np
##################################
# INPUT PARAMETERS:
##################################
sat_file = 'starlette'
arc1 = '030914_2wk'
grav_id ='goco05s'
local_path = '/data/analysis/starlette_analysis/'
SAT_ID = 7501001
##################################
# PATH TO DENSITY MODEL RUN of Choice:
##################################
msis86_model = 'msis86'
path_to_msis86 = '/data/runs_geodyn/st/results/'+ grav_id+'_'+msis86_model+'/'
import sys
sys.path.insert(0, '/data/analysis/util_funcs/py_starlette/')
from a_ReadStarlette import ReadStarlette
AdjustedParams, Trajectory, Density, Resids = ReadStarlette(arc1,
sat_file,
grav_id,
local_path,
path_to_msis86)
The base file name for this arc is: st030914_2wk.goco05s
File exists: iieout
/data/runs_geodyn/st/results/goco05s_msis86/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
/data/runs_geodyn/st/results/goco05s_msis86/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
/data/runs_geodyn/st/results/goco05s_msis86/DENSITY/st030914_2wk.goco05s
Loading data...
Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded
Adjusted Parameters¶
The Adjusted Parameters are read in from the
iieoutfile (Unit 6).These include the State Vector at the start of the ARC, drag coefficient (CD), the solar radiation coefficient (CR), and the adjusted empirical accelerations (GA).
the CD is controlled/output by the DRAG card.
the GA is ouput by the ACCEL9 card (the 9 parameter general acceleration)
AdjustedParamsis a dictionary in which the first dimension stores each geodyn run iteration:
[2]:
AdjustedParams.keys()
[2]:
dict_keys([1, 2, 3, 4, 5])
And the 2nd dimension stores all of the adjusted parameters:
[3]:
AdjustedParams[5].keys()
[3]:
dict_keys(['0CD', '0GA', '0GA_t2', '0XPOS', '0YPOS', '0ZPOS', '0XVEL', '0YVEL', '0ZVEL'])
Each adjusted parameter then has the associated statistics:
APRIORI VALUE
PREVIOUS VALUE
CURRENT VALUE
TOTAL DELTA
CURRENT DELTA
APRIORI SIGMA
CURRENT SIGMA
[4]:
AdjustedParams[5]['0XPOS']['CURRENT_VALUE']
[4]:
-1228610.639773413
In the case of CD and GA, there are multiple values for each, and an additional dictionary exists for each value
[5]:
print(AdjustedParams[5]['0GA'].keys(),'\n')
print(np.squeeze(list(AdjustedParams[5]['0CD'].keys())))
dict_keys(['0GA 9P 11', '0GA 9P 12', '0GA 9P 21', '0GA 9P 22'])
[Timestamp('2003-09-14 08:00:00') Timestamp('2003-09-14 16:00:00')
Timestamp('2003-09-15 00:00:00') Timestamp('2003-09-15 08:00:00')
Timestamp('2003-09-15 16:00:00') Timestamp('2003-09-16 00:00:00')
Timestamp('2003-09-16 08:00:00') Timestamp('2003-09-16 16:00:00')
Timestamp('2003-09-17 00:00:00') Timestamp('2003-09-17 08:00:00')
Timestamp('2003-09-17 16:00:00') Timestamp('2003-09-18 00:00:00')
Timestamp('2003-09-18 08:00:00') Timestamp('2003-09-18 16:00:00')
Timestamp('2003-09-19 00:00:00') Timestamp('2003-09-19 08:00:00')
Timestamp('2003-09-19 16:00:00') Timestamp('2003-09-20 00:00:00')
Timestamp('2003-09-20 08:00:00') Timestamp('2003-09-20 16:00:00')
Timestamp('2003-09-21 00:00:00') Timestamp('2003-09-21 08:00:00')
Timestamp('2003-09-21 16:00:00') Timestamp('2003-09-22 00:00:00')
Timestamp('2003-09-22 08:00:00') Timestamp('2003-09-22 16:00:00')
Timestamp('2003-09-23 00:00:00') Timestamp('2003-09-23 08:00:00')
Timestamp('2003-09-23 16:00:00') Timestamp('2003-09-24 00:00:00')
Timestamp('2003-09-24 08:00:00') Timestamp('2003-09-24 16:00:00')
Timestamp('2003-09-25 00:00:00') Timestamp('2003-09-25 08:00:00')
Timestamp('2003-09-25 16:00:00') Timestamp('2003-09-26 00:00:00')
Timestamp('2003-09-26 08:00:00') Timestamp('2003-09-26 16:00:00')
Timestamp('2003-09-27 00:00:00') Timestamp('2003-09-27 08:00:00')
Timestamp('2003-09-27 16:00:00') Timestamp('2003-09-28 00:00:00')]
[6]:
AdjustedParams[5]['0GA']['0GA 9P 11']
[6]:
{'APRIORI_VALUE': 0.0,
'PREVIOUS_VALUE': 5.536419290656419e-11,
'CURRENT_VALUE': 5.546253541325767e-11,
'TOTAL_DELTA': -5.54625354133e-11,
'CURRENT_DELTA': -9.83425066935e-14,
'APRIORI_SIGMA': 1.0,
'CURRENT_SIGMA': 1.6930522e-11}
[7]:
AdjustedParams.keys()
[7]:
dict_keys([1, 2, 3, 4, 5])
And the 2nd dimension stores all of the adjusted parameters:
[8]:
AdjustedParams[5].keys()
[8]:
dict_keys(['0CD', '0GA', '0GA_t2', '0XPOS', '0YPOS', '0ZPOS', '0XVEL', '0YVEL', '0ZVEL'])
[9]:
import plotly.graph_objects as go
import numpy as np
Cd_plot = []
labels = list(AdjustedParams[5]['0CD'].keys())
for i in AdjustedParams.keys():
Cd_plot.append(float(AdjustedParams[i]['0CD'][labels[0]]['CURRENT_VALUE']))
fig = go.Figure(data=[go.Scatter(x=list(AdjustedParams.keys()),
y=Cd_plot,
mode='markers',
marker=dict(
size=10,
),
),
],
)
fig.update_layout(
title="Drag Coefficient vs Iteration",
yaxis_title=r"$C_d \,\,\text{(Unitless)}$",
xaxis_title="Iteration Number",
)
fig.show()
[10]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
from plotly.subplots import make_subplots
last_iter = 5
which_stat = 'CURRENT_VALUE' # 2 is the current val
labels = list(AdjustedParams[5]['0CD'].keys())
fig = make_subplots(
rows=1, cols=1,
subplot_titles=("Time Dependent Drag Coefficient",
))
val_list = []
for i in AdjustedParams[last_iter]['0CD'].keys():
val_list.append(AdjustedParams[last_iter]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
y=val_list,
name= 'MSIS 86',
mode='markers',
marker=dict(
size=10,),
), row=1, col=1,
)
fig.update_yaxes( title=r"$C_D (Unitless)$",exponentformat= 'power',row=1, col=1)
fig.update_xaxes( title="Iteration Number", row=1, col=1)
# fig.update_yaxes( yaxis_title=r"$\frac{m}{s^2}$", xaxis_title="Iteration Number",exponentformat= 'power',row=2, col=1)
# fig.update_layout(
# title="Along Track Accleration Amplitude vs Iteration", )
# fig.update_layout(
# yaxis_title=r"$\frac{m}{s^2}$",
# xaxis_title="Iteration Number" )
iplot(fig)
Trajectory Data¶
The trajectory data can come from multiple sources. The official output should come from the binary orbit file, but that is not working as of January 2020. Instead, we use the orbit file printout (ORBTVU card) to printout the trajectory data to an ascii file (
Unit 10for keplerian andUnit 08for cartesian). These are then placed in theKEP_TRAJandXYZ_TRAJdata directories before the temporary files are purged.The
Trajectorystructure is a dictionary in which the first dimension represents the different satellites, and the second dimension is a Pandas DataFrame containing the trajectory data from theORBTVUcard in which the first dimension stores each geodyn run iteration:
Indexing the first dimension will show the satellites in the GEODYN run that were saved. They are saved according to their SAT_ID. Since the Starlette-SLR runs only contain the Starlette satellite, this is the only one saved.
[11]:
Trajectory.keys()
[11]:
dict_keys([7501001])
The trajectory data can then be indexed by SAT_ID to show the pandas dataframe containing the data. The stored variables can be found with the .columns method.
[12]:
print(Trajectory[7501001].columns)
Index(['YYMMDD', 'HHMM', 'SECONDS', 'X', 'Y', 'Z', 'XDOT', 'YDOT', 'ZDOT',
'LAT', 'LONG', 'HEIGHT', 'timeHHMM', 'year', 'month', 'day', 'hours',
'minutes', 'secs', 'Date'],
dtype='object')
[13]:
Trajectory[7501001].head()
[13]:
| YYMMDD | HHMM | SECONDS | X | Y | Z | XDOT | YDOT | ZDOT | LAT | LONG | HEIGHT | timeHHMM | year | month | day | hours | minutes | secs | Date | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2 | 30914 | 105 | 55.205719 | -3902733.45 | 5158963.35 | -3130061.00 | -5675.94599 | -1433.45008 | 4725.70967 | -25.954786 | 118.024110 | 812276.80 | 0105 | 2003 | 09 | 14 | 01 | 05 | 55 | 2003-09-14 01:05:55 |
| 3 | 30914 | 106 | 15.005953 | -4014288.94 | 5129496.49 | -3035836.70 | -5591.75463 | -1542.86306 | 4791.45070 | -25.120244 | 118.880387 | 811968.94 | 0106 | 2003 | 09 | 14 | 01 | 06 | 15 | 2003-09-14 01:06:15 |
| 4 | 30914 | 106 | 47.006428 | -4190987.23 | 5077312.08 | -2880862.22 | -5450.70994 | -1718.31837 | 4893.42231 | -23.758977 | 120.237767 | 811607.96 | 0106 | 2003 | 09 | 14 | 01 | 06 | 47 | 2003-09-14 01:06:47 |
| 5 | 30914 | 107 | 13.806875 | -4335433.07 | 5029307.46 | -2748618.90 | -5327.95616 | -1863.82658 | 4974.68198 | -22.607839 | 121.350810 | 811436.43 | 0107 | 2003 | 09 | 14 | 01 | 07 | 13 | 2003-09-14 01:07:13 |
| 6 | 30914 | 107 | 47.207484 | -4510756.72 | 4964053.31 | -2580843.47 | -5169.22701 | -2043.14535 | 5070.56125 | -21.160247 | 122.709738 | 811390.75 | 0107 | 2003 | 09 | 14 | 01 | 07 | 47 | 2003-09-14 01:07:47 |
[14]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
# )
x = Trajectory[7501001]['Z'].astype(float)[:]
y = Trajectory[7501001]['X'].astype(float)[:]
z = Trajectory[7501001]['Y'].astype(float)[:]
fig = go.Figure(data=[go.Scatter3d(x=x,
y=y,
z=z,
mode='markers',
marker=dict(
size=2,
opacity=0.7),
),
],
)
fig.update_layout(
title="3D Cartesian Orbit",
xaxis_title="X Coordinate",
yaxis_title="Y Coordinate",
# zaxis_title="z Coordinate",
)
# fig.show()
iplot(fig)
Density Data¶
GEODYN does not actually have a built in method for printing out the density. As such, the density has been printd by modifying the
DRAG.f90subroutine in theIIE/MODSdirectory. The modification consists of printing the density and trajectory data at each timestep of the GEODYN run.The density data is stored as an ascii file named
fort.99, which is then stored to theDENSITYdirectory before the temporary files are purged.The density data saves the longitutde, latitude, altitude coordinates, and cartesian coordinates.
[15]:
Density.columns
[15]:
Index(['Elapsed Secs', 'YYMMDD', 'HHMMSS', 'Lat', 'Lon', 'Height (meters)',
'rho (kg/m**3)', 'delta_rho (kg/m**3/m)', 'X', 'Y', 'Z', 'XDOT', 'YDOT',
'ZDOT', 'timeHHMMSS', 'year', 'month', 'day', 'hours', 'minutes',
'secs', 'Date'],
dtype='object')
[16]:
Density.head()
[16]:
| Elapsed Secs | YYMMDD | HHMMSS | Lat | Lon | Height (meters) | rho (kg/m**3) | delta_rho (kg/m**3/m) | X | Y | ... | YDOT | ZDOT | timeHHMMSS | year | month | day | hours | minutes | secs | Date | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 210.0 | 30914 | 434 | 45.6973 | 278.5740 | 1083058.594 | 1.623031e-15 | 0.0 | 182992.974 | -5215697.130 | ... | -1841.620957 | -1912.851416 | 000434 | 2003.0 | 9.0 | 14.0 | 0.0 | 4.0 | 34.0 | 2003-09-14 00:04:34 |
| 1 | 225.0 | 30914 | 449 | 45.3666 | 279.6052 | 1084271.860 | 1.606673e-15 | 0.0 | 284238.712 | -5242756.035 | ... | -1766.174236 | -1989.534319 | 000449 | 2003.0 | 9.0 | 14.0 | 0.0 | 4.0 | 49.0 | 2003-09-14 00:04:49 |
| 2 | 240.0 | 30914 | 504 | 45.0256 | 280.6232 | 1085455.504 | 1.590705e-15 | 0.0 | 385422.960 | -5268680.634 | ... | -1690.382976 | -2065.747643 | 000504 | 2003.0 | 9.0 | 14.0 | 0.0 | 5.0 | 4.0 | 2003-09-14 00:05:04 |
| 3 | 255.0 | 30914 | 519 | 44.6746 | 281.6279 | 1086609.356 | 1.575121e-15 | 0.0 | 486523.867 | -5293465.880 | ... | -1614.263040 | -2141.476196 | 000519 | 2003.0 | 9.0 | 14.0 | 0.0 | 5.0 | 19.0 | 2003-09-14 00:05:19 |
| 4 | 270.0 | 30914 | 534 | 44.3139 | 282.6193 | 1087733.252 | 1.559916e-15 | 0.0 | 587519.621 | -5317106.961 | ... | -1537.830323 | -2216.704901 | 000534 | 2003.0 | 9.0 | 14.0 | 0.0 | 5.0 | 34.0 | 2003-09-14 00:05:34 |
5 rows × 22 columns
[17]:
data_nums = 1000
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
fig = go.Figure(data=[go.Scatter(x=Density['Date'][:data_nums],
y=Density['rho (kg/m**3)'][:data_nums],
mode='markers',
marker=dict(
size=2,
),
),
],
)
fig.update_layout(
title="Density along Starlette Orbit (for 4 minutes of first day)",
yaxis_title=r'$\frac{kg}{m^3}$',
xaxis_title="Date",
)
fig.update_yaxes(type="log", exponentformat= 'power',)
iplot(fig)
[18]:
data_nums = 10
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
fig = go.Figure(data=[go.Scatter(x=Density['Date'][::data_nums],
y=Density['rho (kg/m**3)'][::data_nums],
mode='markers',
marker=dict(
size=2,
),
),
],
)
fig.update_layout(
title="Density along Starlette Orbit (Every 10th datapoint)",
yaxis_title=r'$\frac{kg}{m^3}$',
xaxis_title="Date",
)
fig.update_yaxes(type="log", exponentformat= 'power',)
iplot(fig)
Residual Data¶
Similar to the trajectory data, the official output of the residuals should come from the binary residual file (Unit 19), but that is not working as of January 2020. Instead, we read in the residuals that are printed in ascii formats to the
iieoutfile using theOBSVUcard.The residuals are saved as a Pandas DataFrame.
[19]:
Resids
[19]:
| index | YYMMDD | HHMM | Sec-UTC-R | Observation | Residual | Ratio to sigma | Elev1 | Elev2 | OBS No. | ... | timeHHMM | year | month | day | hours | minutes | secs | microsecs | Date | Ratio_to_sigma_fixed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 30914 | 105 | 55.208831 | 933.052599576 | -0.001484 | -0.0015 | 58.485935 | 58.485839 | 1 | ... | 0105 | 2003 | 09 | 14 | 01 | 05 | 55 | 208831 | 2003-09-14 01:05:55.208831 | -0.0015 |
| 1 | 1 | 30914 | 106 | 15.009306 | 1005.251863666 | -0.000437 | -0.0004 | 51.260530 | 51.260437 | 2 | ... | 0106 | 2003 | 09 | 14 | 01 | 06 | 15 | 009306 | 2003-09-14 01:06:15.009306 | -0.0004 |
| 2 | 2 | 30914 | 106 | 47.010250 | 1145.739853454 | -0.002657 | -0.0027 | 41.550447 | 41.550361 | 3 | ... | 0106 | 2003 | 09 | 14 | 01 | 06 | 47 | 010250 | 2003-09-14 01:06:47.010250 | -0.0027 |
| 3 | 3 | 30914 | 107 | 13.811144 | 1279.798448021 | -0.006932 | -0.0069 | 35.007633 | 35.007553 | 4 | ... | 0107 | 2003 | 09 | 14 | 01 | 07 | 13 | 811144 | 2003-09-14 01:07:13.811144 | -0.0069 |
| 4 | 4 | 30914 | 107 | 47.212357 | 1461.073729804 | -0.011274 | -0.0113 | 28.408244 | 28.408171 | 5 | ... | 0107 | 2003 | 09 | 14 | 01 | 07 | 47 | 212357 | 2003-09-14 01:07:47.212357 | -0.0113 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 3999 | 5688 | 30927 | 2306 | 16.259539 | 1260.651605954 | 0.006996 | 0.0070 | 53.959755 | 53.959671 | 4000 | ... | 2306 | 2003 | 09 | 27 | 23 | 06 | 16 | 259539 | 2003-09-27 23:06:16.259539 | 0.0070 |
| 4000 | 5689 | 30927 | 2306 | 46.014220 | 1374.565632664 | 0.002611 | 0.0026 | 46.639968 | 46.639884 | 4001 | ... | 2306 | 2003 | 09 | 27 | 23 | 06 | 46 | 014220 | 2003-09-27 23:06:46.014220 | 0.0026 |
| 4001 | 5690 | 30927 | 2307 | 14.118743 | 1496.430957455 | -0.001994 | -0.0020 | 40.669100 | 40.669018 | 4002 | ... | 2307 | 2003 | 09 | 27 | 23 | 07 | 14 | 118743 | 2003-09-27 23:07:14.118743 | -0.0020 |
| 4002 | 5691 | 30927 | 2307 | 46.624060 | 1650.422940823 | -0.006786 | -0.0068 | 34.755223 | 34.755145 | 4003 | ... | 2307 | 2003 | 09 | 27 | 23 | 07 | 46 | 624060 | 2003-09-27 23:07:46.624060 | -0.0068 |
| 4003 | 5692 | 30927 | 2308 | 5.727222 | 1745.985693706 | -0.003580 | -0.0036 | 31.700982 | 31.700906 | 4004 | ... | 2308 | 2003 | 09 | 27 | 23 | 08 | 5. | 27222 | 2003-09-27 23:08:05.027222 | -0.0036 |
4004 rows × 21 columns
[20]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
fig = go.Figure(data=[go.Scatter(x=pd.to_datetime(Resids['Date']),
y=Resids['Residual'].values.astype(float),
mode='markers',
marker=dict(
size=4,
),
),
],
)
fig.update_layout(
title="Observation Residuals from Final Iteration",
yaxis_title="Residuals",
xaxis_title="Date",
)
iplot(fig)